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One of the most exciting potential sources of gravitational waves for low-frequency, space-based 
gravitational wave (GW) detectors such as the proposed Laser Interferometer Space Antenna (LISA) 
is the inspiral of compact objects into massive black holes in the centers of galaxies. The detection 
of waves from such "extreme mass ratio inspiral" systems (EMRIs) and extraction of information 
from those waves require template waveforms. The systems' extreme mass ratio means that their 
waveforms can be determined accurately using black hole perturbation theory. Such calculations are 
computationally very expensive. There is a pressing need for families of approximate waveforms that 
may be generated cheaply and quickly but which still capture the main features of true waveforms. 
In this paper, we introduce a family of such "kludge" waveforms and describe ways to generate 
them. Different kinds of "kludges" have already been used to scope out data analysis issues for 
LISA. The models we study here are based on computing a particle's inspiral trajectory in Boyer- 
Lindquist coordinates, and subsequent identification of these coordinates with flat space spherical 
polar coordinates. A gravitational waveform may then be computed from the multipole moments 
of the trajectory in these coordinates, using well known solutions of the linearised gravitational 
perturbation equations in flat spacetime. We compute waveforms using a standard slow-motion 
quadrupole formula, a quadrupole/octupole formula, and a fast-motion, weak-field formula originally 
developed by Press. We assess these approximations by comparing to accurate waveforms obtained 
by solving the Teukolsky equation in the adiabatic limit (neglecting GW backreaction). We find 
that the kludge waveforms do extremely well at approximating the true gravitational waveform, 
having overlaps with the Teukolsky waveforms of 95% or higher over most of the parameter space 
for which comparisons can currently be made. Indeed, we flnd these kludges to be of such high 
quality (despite their ease of calculation) that it is possible they may play some role in the flnal 
search of LISA data for EMRIs. 



I. INTRODUCTION-MOTIVATION 



The proposed Laser Interferometer Space Antenna (LISA) j|, is expected to provide a variety of high-precision 
gravitational wave (GW) measurements. One of the most interesting targets for this space-based detector are the 
GWs generated by stellar- mass compact objects inspiralling into (super)massive black holes [(S)MBHs^]. Accumulated 
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^ LISA'S sensitivity is primarily to events involving black holes on the low end of the mass spectrum seen in many galaxies — around 10^ 
to (a few)xlO'^ Mq. We emphasize this point because the name "supermassive black hole" is often taken to refer to a hole with mass 
in the range 10^ — 10® Mq. GW events from such black holes will typically be at frequencies too low for LISA's sensitivity. 
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astrometric observations provide strong support in favor of the existence of a "dark" compact object in the core of 
every galaxy (for which the central parsec region can be resolved) 0. With masses ranging between 10^ — 10® Mq 
these objects are believed to be massive Kerr black holes Q. It is also believed that these holes are the "quiet" 
remnants of an older quasar population ^^l. 

Multi-body interactions in the "cusp" stellar population surrounding these SMBHs can put stellar mass compact 
objects onto orbits that come close to the central black hole. If the object passes sufficiently close to the (S)MBH, 
it may be captured and subsequently spiral in by the emission of GWs [1, 0, 0| • Initially, the captured bodies are 
expected to be on "generic" orbits, i.e., eccentric (with eccentricity e « 1) and inclined with respect to the central 
black hole's equatorial plane 0, Q. These orbits evolve adiabatically due to GW emission, decreasing in eccentricity 
and periastron while the inclination of the orbit remains approximately constant, but increases slightly (see [9l.[l0l[Tl| 
for an approximate description of the full inspiral). 

For central black hole masses in the range IO^Mq < M < lO^M©, the GWs emitted during the inspiral will be at 
frequencies close to the floor of the LISA noise curve (~ 3 mHz). LISA will detect the bursts of radiation emitted 
near periapse throughout the inspiral, but these will not be individually resolvable unless they occur in our 
Galaxy [13]. However, during the last few years of inspiral, when the small object is orbiting deep inside the SMBH's 
gravitational field, GWs are radiated continuously in the LISA band. During this phase, the residual eccentricity is 
typically at e ^ 0.4 [1, [l3| and the orbital motion will exhibit extreme versions of relativistic effects (e.g., periastron 
and Lcnsc-Thirring precession). As a consequence, the resulting GW signal will be strongly "colored" by these effects 
and will take a complicated form [14, 15j. By complicating the waveform, these strong-field effects potentially make 
data analysis difficult; but, they also encode a great deal of information about the strong field nature of the spacetime. 
By accurately measuring all of these effects, it is expected that we will be able to "map" the spacetime of the black 
hole [III , probing its multipolar structure and verifying that it obeys the constraints that general relativity imposes 
on black hole solutions [17|, ll8| . 

The expected amplitude of the signal from these extreme-mass-ratio inspirals (EMRIs) is about an order of mag- 
nitude below LISA'S projected instrumental noise and, at low frequencies (< 2 mHz), is several orders of magnitude 
below "confusion noise" produced by unresolved Galactic binaries [Toj. However, the signals will be observed for 
~ 10^ waveform cycles, and matched filtering will therefore allow detection of these signals with high signal-to-noise 
ratio (SNR) out to a redshift z ~ 1 [1^. Preliminary estimates suggest LISA could see as many as ~ 10'^ EMRI 
events during its lifetime 0, , using a suitable "semi-coherent" search algorithm and provided that confusion noise 
can be efficiently reduced in the real data. 

Matched filtering algorithms require the correlation of the detector's output data stream with a bank of waveform 
templates which describe the real signal with sufficient accuracy, covering the whole parameter space. The fact that 
we are dealing with a binary system of extreme mass ratio ^/M <C 1 means that the gravitational waveform may be 
obtained accurately using black hole perturbation theory. The extreme mass ratio also guarantees that the orbital 
parameters evolve adiabatically under radiation reaction; i.e., they evolve on a much longer time scale than the 
orbital periods. This implies that within the radiation reaction time scale, the inspiral waveform can be approximated 
by "snapshot" waveforms — waveforms calculated by assuming that the small object is moving along a geodesic, 
neg lecting backreaction for that short stretch of time. These snapshots are constructed using the Teukolsky equation 
[2l| . an equation that describes the first order change to the curvature tensor of a black hole due to some perturbing 



source. Accurate Teukolsky-bascd (TB) snapshot waveforms have been calculated for inclined-circular orbits 14 1, 
equatorial-eccentric orbits [15,, 2M ^i^'i most recently, for a certain number of generic (inclined-eccentric) orbits 124 1 . 
The reader can find recent reviews on the modeling of EMRI waveforms and orbital dynamics in Refs. pTl[25|. 

TB waveforms are computationally expensive to generate as they require the numerical integration of the Teukolsky 
equation and summation over a large number of multipole modes. In this sense, these waveforms are not "user- 
friendly" , especially when it comes to realistic data analysis calculations where one has to handle a bank which 
contains ~ 10"'^^ of these waveforms i2(| . Moreover, the Teukolsky formalism does not provide any information on 
"conservative" self-interaction effects. To compute these, one must use a self-force formalism. This approach is still 
under development and is very likely to be even more computationally expensive when it is completed (see p6j for a 
recent review). 

These difficulties have led to a need for the construction of approximate families of waveforms that capture the 
main features of the true signals, but which can be generated quickly in large numbers. Such approximate waveforms 
are already being used for scoping out data analysis algorithms for the detection of EMRIs in LISA data [20j , and 
may ultimately play some role as fiducial detection templates in the final data analysis. One possible approach is 
to construct post-Newtonian (PN) waveforms, which have the advantage of being analytic and therefore very easy 
to generate. Post-Newtonian EMRI waveforms have been computed in the Schwarzschild spacetime for both quasi- 
circular [2711 and eccentric orbits [28] , and in the spacetime outside a slowly rotating Kerr black hole for quasi-circular 
orbits [231 . 

Recently another class of approximate waveforms have become available, based on various "kludge" approaches [sot . 
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The basic idea of the kludges is to combine different prescriptions for the orbital evolution and waveform emission 
(not necessarily in a self-consistent way). An "analytic kludge" (AK) was developed by Barack and Cutler [3l| (see 
also [i^])- In their model, the small object is moving on a Keplerian orbit, amended to include the effects of pericenter 
precession, Lense-Thirring precession, and inspiral from radiation reaction. The emitted gravitational waveforms are 
described by the lowest-order quadrupole formula. 

In this paper, we consider an alternative way to construct kludge waveforms. This approach is much less amenable 
to analytic calculation, so we called it the "numerical kludge" (NK). The first step of the NK is to produce an inspiral 
trajectory in "phase space" — the space defined by the constants (orbital energy, axial angular momentum, and 
"Carter constant" ) which specify orbits (up to initial conditions) d, . The second step is to numerically integrate 
the Kerr geodesic equations along this insp iral trajectory and hence obtain the Boyer-Lindquist coordinates of the 
inspiralling object as a function of time |33| . The final step is to construct a gravitational waveform from this inspiral 
trajectory. 

The approach to waveform construction that we take is to identify the Boyer-Lindquist coordinates of the source 
with spherical polar coordinates in flat-space. There are several different expressions available in the literature for the 
gravitational waveforms from orbits in flat space, and we apply these to our pseudo-flat-space trajectory to construct 
waveforms. Specifically, we look at the quadrupole formula ^34] , which is valid in the limit of weak-field (i.e., flat- 
space) and slow motion. We also examine the quadrupole-octupole formula of Bekenstein [3^ , as well as a formula 
derived by Press [1^ . The Press formula is also a weak-field expression, but is not restricted to slow- motion or small 
sources, and contains radiation at orders higher than quadrupole and octupole. 

The purpose of this paper to establish and delimit the accuracy and reliability of waveforms constructed in these 
various ways. We do this by comparing to TB waveforms. TB waveforms are currently the most accurate EMRI 
waveforms available. In most cases, TB waveforms represent the emission from geodesic orbits — we mostly do 
not include the radiative evolution of orbital parameters in this work. (With one exception: because complete TB 
inspirals exist for zero-eccentricity orbits j37| . we compare to a full inspiral in this case.) We compare the various 
NK waveforms with TB waveforms using an overlap integral which weights the waveforms in frequency space by the 
expected LISA noise curve, and maximizes the overlap with respect to time offsets. This overlap is identical to the 
test used to evaluate the efficiency of model waveforms as detection templates. 

We find that quadrupole-octupole NK waveforms are in remarkable agreement (overlaps ^ 0.95) with TB waveforms 
for orbits with periastron Vp > 5M . Most orbits in the final year of the inspiral satisfy this restriction which means that 
NK waveforms are quite accurate over a considerable portion of the inspiral parameter space. For orbits that come 
below that radius (mostly prograde orbits around rapidly spinning holes), the agreement rapidly degrades (although 
NK waveforms remain more accurate than post-Newtonian or AK waveforms). This is not surprising — for such 
orbits the full TB waveform receives significant contributions from higher multipoles and back-scattering from the 
spacetime. These effects are ab initio absent in the NK prescription. Waveforms generated using the Press formula 
perform similarly well for Tp ^ 5M, and slightly better than the quadrupole-octupole prescription for orbits where the 
contribution from higher harmonics is not negligible. The improvement in this relatively strong-field regime is only 
slight, however — we did not find any regime where the Press formula is a significant improvement on the quadrupole- 
octupole prescription. We conclude that NK waveforms can accurately reproduce true gravitational waveforms in a 
large part of parameter space; but, there is little gain from going beyond the quadrupole-octupole prescription. 

We also calculate kludge radiation fluxes by combining the quadrupole energy flux formula with exact geodesic 
motion and averaging in time; such an analysis was done for Schwarzschild orbits in ^38J . The resulting fluxes 
typically compare well to the best available PN formula in the weak field, but provide a somewhat better estimate 
of the energy and angular momentum fluxes (comparing to TB results) for strong-field orbits. These fluxes also 
allows us to assess by how much the approximations which go into the kludge construction are inconsistent. Since 
the inspiral trajectory and waveform construction are considered separately, kludge GWs carry a different amount 
of energy and angular momentum to infinity than the inspiralling particle loses (which is set by the formulas which 
determine the inspiral through orbital parameter space). It is important to bear this in mind when using kludge 
waveforms in computations, e.g., for estimating SNRs of LISA EMRI detections, and to have an estimate of the 
size of the inconsistency. Finally, we describe how the kludge waveforms may be improved by including some of the 
conservative self-force effects, which may have a significant imprint on true inspiral waveforms. 

This paper is organised as follows. Section |TT] provides a review of existing EMRI waveforms, paying special 
attention to the kludge semi-relativistic waveforms. Details of how these waveforms can be generated are described. 
In Section Hill we discuss the overlap function between two waveforms, as this will be our main tool for comparison. 
Section llVI contains the results from the comparison between kludge and Teukolsky-based waveforms and fluxes. In 
Section |V] a method is outlined for including the conservative parts of the self- force in our kludge scheme. Finally, we 
present a concluding discussion in Section IVTl We shall use greek letters (/x, i/, etc.) to denote spacetime indices, and 
latin letters (i, j, etc.) to denote spatial indices (unless explicitly stated otherwise). Throughout the paper we adopt 
geometric units G = c ~ 1. 
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II. WAVEFORM INVENTORY 



Presently, several types of EMRI waveforms are available. In broad terms, these waveforms fall into three categories: 
(i) those calculated numerically within the framework of black hole perturbation theory (Teukolsky-based, or TB 
waveforms) , (ii) analytic waveforms which result from self-consistent PN expansions of the GW equations and orbital 
motion, and (iii) approximate semi-relativistic waveforms, or "kludges". This third category is the focus of this paper. 

Kludge waveforms are constructed by combining a flat spacetime wave-emission formula with either a fully relativis- 
tic treatment of particle motion (resulting in the numerical kludge, or NK, waveforms), or some analytic approximation 
of this motion (leading to the analytic kludge, or AK, waveforms). We shall examine the construction of NK wave- 
forms, comparing them to TB waveforms (the most accurate EMRI waveforms presently available). For completeness 
and background to this paper's discussions, we now briefly discuss each of the above waveform families. 



The primary framework for black hole perturbation theory in a Kerr background is the Teukolsky formalism [2l|, 
which encapsulates all gravitational radiative degrees of freedom in a single "master" wave equation — the "Teukolsky 
equation" — for the Weyl scalars "00 and "04. A key feature of this equation is that it admits separation of variables 
in the frequency domain, which effectively reduces it to a pair of ordinary differential equations. There are extensive 
results in the literature on solutions of the Teukolsky equation, starting from Press and Teukolsky [s^ , Detweiler 40 1 
and Sasaki and Nakamura [4l|. More recent work uses the Teukolsky-Sasaki-Nakamura formalism, see Refs. [ll|, 25 1 
for detailed discussions and references on the subject. 

To date, the Teukolsky equation has been solved in the frequency domain for circular-inclined orbits p^ , eccentric- 
equatorial orbits [3] and most recently for a handful of generic (eccentric-inclined) orbits [24]. We make use of the 
waveforms generated by these various authors in this paper to assess the quality of our NK waveforms. Recently, the 
Teukolsky equation has also been solved directly in the time domain [Z^, Sj] . Time-domain calculations have the great 
advantage of speed, since they avoid the need for Fourier decomposition and summation over all orbital frequency 
harmonics. However, these calculations are not yet mature enough to provide accurate waveforms from Kerr orbits, 
because of the difficulty of representing the various (5-functions appearing in the energy-momentum tensor of a point 
particle. 

As a final remark, we should mention again that in all the above TB calculations (either in the frequency or the 
time domain) the motion of the small object is taken to be strictly geodesic. This is a reasonable first approximation 
since for an EMRI system the orbital evolution is adiabatic, i.e., over a time interval of several orbits the motion is 
geodesic to high precision. However, for longer stretches of time (~ Af^//i) the effects of gravitational back-reaction 
become significant and cannot be ignored anymore. Waveforms that take into account an evolving orbit (and the 
conservative self-interaction) require the computation of the gravitational self-force (see [43 | for an up-to-date review 
and a full list of references). However, self- force waveforms are not yet available and are unlikely to be for the next 
few years. Moreover, it is very likely that self-force calculations will remain quite computationally intensive; as such, 
it is very likely that it will not be possible to generate self-force based waveforms in sufficient numbers to be used for 
LISA data analysis. It is therefore essential to investigate approximate, easy-to-use waveform models. 



Most available analytic waveforms are based on post-Newtonian expansions of the orbital dynamics and wave 
emission, an expansion that is of greatest relevance when the bodies are widely separated. These waveforms are 
typically constructed for a specific object or for restricted orbital motion. A significant amount of work has been 
done on modeling waveforms from two spinning bodies with comparable masses orbiting in precessing quasi-circular 
orbits (a key GW source for both ground- and space-based detectors). Kidder [451 has investigated the effects of 
spin-spin and spin-orbital coupling on the waveform from inspiralling binaries. The most recent investigation of PN 
waveforms (and their application to data analysis) for spinning binary systems may be found in Refs. [i^, 113, SI] 
(see also references therein). Another promising approach is described in [i^ . which is the first attempt to extend 
the "effective-one-body" method [50| to spinning compact objects. In parallel to modeling spinning binaries in quasi- 
circular orbits, there has been significant progress in the construction of waveforms for eccentric, comparable mass 
binaries [sil . [53 . [53| . Presently, there are no post-Newtonian waveforms which include both eccentricity of the orbit 
and spins of the orbiting bodies. 

Post-Newtonian models are ultimately not so useful for modeling EMRIs, since most of the GWs observable to 
LISA are generated from a strong field region (r ^ a few M), where the PN expansion is unlikely to be reliable. One 



A. Teukolsky-based numerical waveforms 




B. Analytic waveforms 



5 



can, however, construct PN waveforms in the EMRI hmit, accepting their certain unreHabihty as a way to develop 
a very "quick and dirty" set of tools for studying these waves. PN EMRI waveforms are available for systems of a 
test mass in a quasi-circular [131 or eccentric [13] orbit around a Schwarzschild black hole, or in a quasi-circular orbit 
around a slowly rotating Kerr black hole [2§| . 

More recently, a class of approximate PN waveforms has been developed by Barack and Cutler [31]. These "ana- 
lytic kludge" (AK) waveforms are essentially phenomenological waveforms — they are constructed using the classic 
quadrupole waveforms for eccentric Keplerian orbits derived in [5^, but with the relativistic effects of pericenter 
precession, Lense-Thirring precession, and inspiral imposed. Though not as accurate as the NK waveforms described 
in this paper, the AK are very quick to generate, and have proven to be useful for computing the Fisher information 
matrix in investigations of parameter measurement with EMRI GWs [3l|. 

The overlap between AK and NK waveforms is high in the very weak field, but degrades as the orbital periapse 
is decreased (see (20| : more detailed comparisons will be included in future papers on the semi-coherent algorithm 
currently in preparation). Even for geodesic orbits, AK and NK waveforms with the same physical parameters drift 
out of phase quickly, since the frequency structure of the two waveform families is different. This arises because the 
AK uses a Keplerian orbital parameterisation, compared to the geodesic parameterisation used in the NK. For an 
equatorial orbit with semi-latus rectum p = 30M and eccentricity e = 0.3 around a IO^A/q black hole of spin a = 0.8, 
the azimuthal frequency of the NK waveform is 0.196mHz compared to 0.216mIIz for the AK. These orbits will 
therefore be one cycle out of phase within ^ 6 hours. AK waveforms will thus not be particularly faithful templates. 
The problems can be mitigated by adjusting the orbital parameters of the AK waveform to improve the match with 
the NK, and the AK waveforms do capture the main features of true EMRI waveforms. For this reason, they may 
be quite effectual templates, but this has not yet been properly assessed. In the future, the effectualness of the AK 
waveforms as detection templates will be investigated by using banks of AK templates to search for more accurate 
NK or TB waveforms embedded in noise. 



C. Semi-relativistic numerical "kludge" waveforms 

The idea behind the numerical kludge (NK) waveforms is to combine an exact particle trajectory (up to inaccuracies 
in the phase space trajectory and conservative radiation reaction terms) with an approximate expression for the GW 
emission. By including the particle dynamics accurately, we hope to capture the main features of the waveform 
accurately, even if we are using an approximation for the waveform construction. 

The computation of NK waveforms proceeds in three steps. The first is to construct an inspiral trajectory in 
"phase space" — that is, the space of constants E (energy), (axial angular momentum) and Q (Carter constant) 
which characterizes Kerr black hole orbits (up to initial conditions). The construction of the phase space trajectory 
has already been described in previous work [§, [l3| ■ In this paper we shall (in most cases) ignore the evolution of 
orbits due to radiation reactions. The procedure for waveform construction including orbital evolution is identical 
to that given here. The second step is to integrate the Kerr geodesic equations along the inspiral trajectory and 
hence obtain the Boyer-Lindquist coordinates of the inspiralling object as a function of time |33j . The final step is 
to construct a gravitational waveform from this inspiral trajectory. We do this by identifying the Boyer-Lindquist 
coordinates (r, 9, 0, t) with spherical polar coordinates in flat-space and then evaluating a flat-space emission formula 
for the corresponding flat-space source orbit. 

The idea of coupling a weak field formula with fully relativistic motion first appeared in papers by Sasaki and 
Ruffini [55I |56| and was then investigated more thoroughly for a test-body orbiting a Schwarzschild black hole by 
Tanaka et al. [23|. More recently, Gair et al. 38] looked at semi-relativistic fluxes for Schwarzschild orbits and derived 
analytic expressions for the fluxes from arbitrary orbits. Their focus was on highly eccentric orbits relevant to the 
capture problem. In all these cases, the authors were interested in computing semi-relativistic fluxes, rather than 
waveforms. While our focus is on waveforms rather than fluxes, in Section FlV CI we do calculate semi-relativistic fluxes 
for Kerr orbits, as these results have not yet appeared in the literature. 

In the remainder of this section, we describe the two stages of waveform generation (ignoring backreaction): (i) the 
computation of a trajectory, and (ii) the computation of a gravitational waveform from an arbitrary trajectory. 



1. Computing the orbital trajectory 

The first step in constructing a numerical kludge waveform is to compute the trajectory that the inspiralling body 
follows in the Boyer-Lindquist coordinates of the Kerr spacetime of the central black hole. Ignoring radiation reaction, 
this path is a Kerr geodesic. Geodesic motion in the Kerr space-time is well-known [s^, [Eg and is governed by a set 
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of first-order differential equations: 



dr 
^d^ 


= ±V 
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= ±v 




d(t> 


= ^0, 




dt 
^d^ 


= Vt, 





where the various "potentials" are defined by 



Vr = [E{r^ + a^)- L,a] - A [r^ + {L, - aE'f + Q] , 
Ve = Q- 008^61 a^{l-E^) 



sin^ e A 

,2 Q 



sm 



Vt = a [Lz - aE sin^ 



- Lza] , 
'E{r'^ + a^) - LM 



(la) 
(lb) 
(Ic) 
(Id) 

(2a) 
(2b) 

(2c) 

(2d) 



Here, E 



1^ cos^ 9, and A 



2Mr + a? . The constants E, L^, Q are the three first integrals of the motion: 



E is the orbital energy; Lz is the projection of the orbital angular momentum along the black hole's spin axis; and Q 
is known as the "Carter constant" . This last constant is the relativistic generalization of the "third integral" used to 
separate the equations which describe orbits in an axisymmetric gravitational potential (a result which is particularly 
well-known in the literature describing orbits in galactic potentials [s^). In the spherical (i.e., Schwarzschild) limit, 
Q reduces to the square of the orbital angular momentum projected into the equatorial plane; see [9[ and references 
therein for discussion. 

For a given E, and Q, the roots of Vr determine the turning points of the radial motion — the periastron rp, and 
apastron r^. From these, one can define an orbital eccentricity e, and semi-latus rectum p, using the conventional 
Keplerian definitions 

P P 



1-t-e' 

2 Tp 



1 



Ta + Tp 



We also replace the Carter constant by an "inclination angle" , defined by 



tan^ L 



LI' 



(3a) 
(3b) 

(4) 



To aid numerical integration, one can work in terms of two angular variables, "0 a-nd instead of r and 6. We define 
■0 by the equation 

P 



\ + e cos V' ' 

We define x by the equation z = cos^ 6 = z_ cos^ x, where z_ is given by 

/3(z+ - z)(z_ ~ z)^ (3z^ ~ z[Q + Ll + a^{\ - E"^)] + Q, 
with f3 — a^(l — E'^). Expanding the radial potential as 

Vr - (1 - E^) (r, - r) (r - Vp) (r - rg) (r - u), 
we find evolution equations for ip and x of the form 

9 - - 

dip \J\ — E^ [(p — r3(l — e)) — e(p + r3(l — e) cos V')] ^ [(p — r4(l -|- e)) -|- e(p — r4(l + e) cos 0)] = 



dt 

dx 
At 



[7 + a2i;z(x)](l-e2) 



V/3 [z+ - z(x)] 
7 + E z{x) 
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In terms of the variables (j), ip and x, the geodesic equations are well behaved at the turning points of the motion, 
which facilitates numerical integration (as discussed in [ssjl. 

Although we mostly ignore it here, it is easy to include radiation reaction in this prescription. To do so, one first 
computes an inspiral trajectory through phase space by writing 

— = fE{a,AI,m,p,e,L), (10a) 

^ = fLia,M,m,p,e,L), (10b) 

^ = fQia,M,m,p,e,L). (10c) 

The functions /e, /l and /g are derived from well-known post-Newtonian results, augmented by additional correc- 
tions. The leading order part of Je is given in a later section of this paper, Eq. (|46p . (For the full expressions see 
Ref. 10]). Since the parameters a, M and m are constant, and p, e and t are directly related to E, and Q via 
Eqs. ©-([I]), equations (fTO|) can be integrated to give the phase space trajectory. That is, the values of E', and Q 
(or equivalently p, e and l) are given as functions of time. To obtain the trajectory of the inspiralling particle, these 
time dependent expressions can be substituted into the right hand sides of Eqs. px|) . (fid)) . ^ and The resulting 
set of ODEs can then be integrated to give the inspiral trajectory. In the following (with one exception), we will only 
consider waveforms from geodesic trajectories (setting /e = /l = fg = 0), and refer the reader to Refs. [1, [l3| for 
more details on the construction of inspirals. 

Once a trajectory has been obtained in this manner, one constructs an "equivalent" flat-space trajectory by pro- 
jecting the Boyer-Lindquist coordinates onto a fictitious spherical polar coordinate grid, defining the corresponding 
Cartesian coordinate system and pretending that these new coordinates are true fiat-space Cartesian coordinates: 

X — r sin 6* cos (/), (Ha-) 
y = rsin6'sin0, (lib) 
z — rcos9. (11c) 

We use the resulting flat space trajectory as input to a wave generation formula. This is a "bead on a wire" 
prescription — by putting the trajectory in flat space, we remove the background that is causing the curvature of the 
geodesic path and hence we are forcing the particle to move along a curved path like a bead moving on a wire. This 
leads to obvious inconsistencies in the approach — e.g., the non-conservation of the flat-space energy- momentum 
tensor of the particle motion since we are not including the energy-momentum of "the wire" along which the particle 
moves. 



2. Waveform generation 

Having constructed the particle orbit in our pseudo-flat-space, we now apply a flat-space wave generation formula. 
Consider the weak- field situation and write down the spacetime metric as (7^,^ ~ rj^^j + h^^, where "q^i, is the flat metric 
and hfiu are small perturbations. The trace-reversed vaeiiic perturbation is defined as /i''" = /i'"^ — {l/2)ri'^'^h, where 
h = rj^'^hfj^^. Imposing the Lorentz gauge condition = 0, the linearized Einstein field equations can be written as 

nh'^'' = -IGTrT^'', (12) 

in which □ denotes the usual fiat space wave operator and the effective energy-momentum tensor T^"^ satisfies 

T'"',^ = 0. (13) 

Here a comma subscript denotes partial derivative (/.^i = df/dx^). Note that our source conservation equation 
uses a partial rather than a covariant derivative. This is because we would hope, in a self-consistent approach, to 
choose coordinates so that the energy momentum tensor is flat-space conserved. Finally, when observing GWs at 
large distances, we are really only interested in the transverse and traceless parts of the spatial components of h^^] a 
projection of these components will be necessary. 

Taking coordinates centered at the black hole, we denote the position of the observer by (i, x) and the position of 
the particle by (ip,Xp)). The wave equation (|12p has the familiar retarded time solution 

4/l!^(i^fe^dV^ (14, 
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The additional coordinate x' is the integration variable, which goes over every possible space location where the 
effective energy- momentum tensor T^'^{t' , x') is nonzero. If the source motion is only negligibly influenced by gravity, 
then T'^'^ may be taken to equal T^^'^ , the energy-momentum tensor of the matter source. In Ref. [IBl, Press derived a 
formula valid for extended, fast motion sources. This was obtained by substituting Eq. into Eq. (|14p repeatedly. 
The result is 



niUm x-'x 



J t'=t-|x-x'| 



„ d'x', 



(15) 



in which = x • x and n = x/r. In the slow motion limit, the Press formula reduces to the usual quadrupole formula 
(hereafter, an overdot denotes a time-derivative) 



(16) 



where 



P\t') = / x'^a;'^-T00(t',x')d3x' 



(17) 



is the source's mass quadrupole moment. Including the next order terms (the mass octupole and current quadrupole 
moments of the source, denoted M^^^ and S^^^ respectively), we obtain the quadrupole-octupole formula (see Refs. [H, 
|36I] for details). 



(18) 



with 



7Vfyfe(t') 



x"x'^x"'T"°it',x')d^x'. 



(19) 
(20) 



Note that in both Eqs. (fT6|) and (fT8|) . the retarded time is t~r instead of the more complicated expression appearing 
in (lisp . If desired, it is a straightforward (but increasingly tedious) task to include more terms in the slow-motion 
expansion of In the present work we shall make use of Eqs. and the full Press formula p3|) . 

The waveform in the standard transverse-traceless (TT) gauge is given by the TT projection of the above expressions. 
We define an orthonormal spherical coordinate system via 



d 



1 d 



dr' ^ r 98 ' 



e$ = 



1 



d 



r sin Q 



(21) 



The angles {Q, denote the observation point's latitude and azimuth respectively. The waveform in transverse- 
traceless gauge is then given by 





/lee-Zi** 2/1^* 
2/i«* h'^" - /i^e 



(22) 



with 



h^^ = cos^ e [/i^^ cos^ $ + h'^y sin 2$ + sin^ $] -|- /i^^ sin^ 9 - sin 29 cos $ + /i^^ sin $] , 



= cos 9 



i /i"'"' sin 2$ + cos 2$ -f- i hyy sin 2$ 



sin 9 [/i^^^ sin <i> - ft,^^ cos $] 



/i** = h'=''iiin'^-h'=ys\n2^ + hyycos^^. 



(23a) 
(23b) 
(23c) 



The usual "plus" and "cross" waveform polarizations are given by /i®® — /i** and 2/i®* respectively. The expres- 
sions p^ , pB|) and pB]) are valid for a general extended source in flat-space. If we specialize to the case of a point- mass 
/i moving along a trajectory x'^{t), we know that the energy-momentum tensor in flat spacetime is given by 



T'^'^(t',x') 



-^-^5^ {x' - x'^ir)) dr = f^^^^S^ (x' - x',(0) , 



dr dr 



dtp dr dr 



(24) 



9 



where r is the proper time along the trajectory. It is related to the particle's coordinate time by dr = (1 — w^/c^)^/^ dt', 
where = \ dx^ /dip p . On the right hand side of Eq. ^ is a term dt'^/dr = l + 0{v^ /c^). Expressions ^ and ^ 
are slow-motion expansions and at the order of these expansions we should replace dt^/dr by 1 for consistency. The 
Press formula (jlSp is a fast-motion expression, and therefore we do include this term. The presence of a (5-function in 
T^" facilitates the simplification of the various moments in Eqs. and p8)) : 

P'' - fix'^x'^' , (25) 
gHk ^ jjk ^ (26) 
M'^'' = x'^ . (27) 

Here = dx"^ /dt'^. The Press formula similarly simplifies to 

(28) 

tp=t-|x-x;| 

The square bracketed expression is to be evaluated at a time t'^ given imphcitly hy t — t^-l- |x — x'p(tp)|. In a numerical 
implementation, we evaluate the expression as a function of the time t' along the particle's path. In order to obtain 
a time series that is evenly spaced in t and compute the RHS of Eq. (^8)) by finite-differencing, we adjust the spacing 
of the sampling at the particle, Ji^, such that 

5t={l-7laV'')5t'j,. (29) 

With a delta function source ((24|) . the retarded-time solution (fT4l) can be evaluated directly [s^: 

Naively, one could suppose that ([30]) will perform better than pS]) . since it is derived using only one of the two 
(invalid) flat space equations, rather than both. In fact, we find that the retarded integral expression ([50)) performs 
much worse than either (fTB|) . (US]) or (|15p when compared to TB waveforms. The reason appears to be that the 
manipulations which lead to the quadrupole, quadrupole-octupole, and Press formulae ensure that the actual source 
terms — mass motions — are on the right hand side. This is to be contrasted with the retarded integral expression 
(fT4| . which identifies the dominant GW on the left hand side with weak spatial stresses on the right. 

We must emphasize that the NK prescription is clearly inconsistent — we are binding the particle motion to a Kerr 
geodesic while assuming flat spacetime for GW generation and propagation. This is manifested by the fact that the 
energy-momentum tensor ([24)) is not flat space conserved, d^T^^'^ ^ 0. However, the spirit of this calculation is not 
a formal and consistent approximation to EMRI waveforms; it is rather a "phenomenological" approach which takes 
into account those pieces of physics we believe are the most crucial — in particular, the exact Kerr geodesic motion. 
By including the exact source trajectory, we ensure that the spectral components of the kludge waveforms are at the 
correct frequencies, although their relative amplitudes will be inaccurate. As we shall see, this line of thinking is 
validated post facto by the remarkable agreement between the kludge and TB waveforms (see Section ITV]) . 

It is important to underline the physical assumptions that have been made in the derivation of (jlSp . (|16p and (jl8p . 
in order to understand their generic limitations. First of all, the assumed absence of any background gravitational field 
means that our kludge waveforms are unable to capture any features related to back-scattering. This effect is known 
to first appear at 1.5PN [i.e. 0{v'^)] level (see, e.g.. |60|). Such "tails" of waves are particularly prominent in the 
strong-field TB equatorial "zoom- whirl" waveforms |l5| and in the waveforms from plunging or parabolic orbits [6l| . 
In all these cases, the hole's quasinormal mode ringing adorns the emitted signal. 

The slow-motion nature of expressions (fTB)) and (fTS]) suggests that they might be bad models for waveforms generated 
by orbits venturing deep inside the central BH's spacetime. The rich multipole structure of the true waveform from 
such orbits is poorly reproduced, as any slow-motion approximation essentially truncates the sum over multipoles. 
The Press formula (fT5|) includes contributions at all multipoles and so might be expected to handle these contributions 
quite well. However, it turns out that it does not perform that much better than the quadrupole-octupole waveform 
in this regime. While it includes contributions at all multipoles, the lack of background curvature in the waveform 
model means that the amplitude of the higher modes is much lower than expected for true EMRI waveforms. 

Fortunately, these two deficiencies become important for the same class of orbits — those that allow the body to 
approach very close to the black hole. As a rule of the thumb we shall find that NK waveforms are reliable (in terms 
of the overlap discussed in the next section) as long as the closest orbital approach (periastron) is > 5M. 



h^''{t,^) 



2fi d^ 
Vdt2 



4^ 
r 



dt^ 
dr 



dx'p dx'p 



dt; 



dt; 



UaV" 



(30) 
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III. THE OVERLAP BETWEEN WAVEFORMS 



We now give a brief description of the measure used for the quantitative comparison between NK and TB waveforms. 
The main motivation for the computation of accurate EMRI waveforms is to carry out matched filtering for detection 
of the GWs. For the purpose of signal detection and parameter estimation, a bank of templates is constructed 
which covers the desired parameter space with sufficient resolution. The detector output is then filtered through each 
template. The measured strain amplitude in the detector, x{t) = s{t) + n(t), consists of a (possibly present) signal 
s{t), and the detector noise n{t). We define the Fourier representation of these time series as 

/>OC 

i{f)= / x{t)e-'^^'fdt. (31) 



The signal-to-noise ratio (SNR) can be expressed in terms of an inner product defined on the vector space of possible 
signals. Given two vectors (time series), x{tk)-, h{tk), we define the overlap {x\h) by the equation [63.l63t 

where 5/1 (/) is the one-sided noise power spectral density (PSD) (in this case, the noise PSD for the LISA detector) and 
an asterisk denotes complex conjugation. Considering h{tk, A") as a template with parameters A" we can approximate 
the SNR by [6| 

i[hix'^)^ ^ jm^ ^ aeb= , (33) 

^ Viinmn Vim 

The notation (/)„ means to ensemble average the function / over all possible noise realizations n. The PSD in the 
denominator of Eq. p2p serves to suppress those frequency components of the signal at which the detector noise is 
large. 

The main tool which will be used in this paper is not the SNR, but rather the overlap function O. The overlap is 
defined as an inner product between two normalized vectors/signals: 

O = is\h) . (34) 

The normalization is chosen so that (s|s) = {h\h) = 1. The overlap can be regarded as the inner product between two 
unit vectors; it varies within [—1, 1]. The overlap is equal to 1 if the two waveforms are identical, and it equals zero 
if the two waveforms are orthogonal (for example, cosine and sine signals). The overlap is an appropriate measure 
of "goodness-of-fit" since we are interested in knowing how well the NK approximate the behavior of TB EMRI 
waveforms. In this context, it is important to include the noise properties of the detector since it is no problem if a 
template has poor performance at frequencies where detector noise is large. 

For a fair comparison, we should choose the signal (TB waveform) and template to have the same parameters^. In 
practice, the start times of our signals and templates were slightly mismatched, so we allow for maximisation over 
starting time ("time of arrival"). This maximisation is commonly done in GW data analysis, and can be accomplished 
very efficiently in the Fourier transform, since the time offset corresponds to a phase shift in the frequency domain. 
Accordingly, we modify the inner product {x\h) slightly: 

where t corresponds to the time difference. We find that the maximum usually occurs with an time offset close to 
zero (typically a few bins). 

We conclude our discussion of the overlap by describing our approximation to the LISA noise Sh{f)- We have used 
an analytic approximation to the numerically generated sensitivity curve given in Ref. ,65j . The agreement between 



^ Note that here we study not the fitting factors but the faithfulness of the NK waveforms as compared to TB-based ones, i.e., how well 
an approximate waveform with given parameters reproduces the "true" waveform with the same physical parameters. A faithful bank 
of waveforms could be used for parameter estimation, while for detection all that is required is an "effectual" template bank, i.e., one 
in which every true waveform is well represented by one template, even if that template has very different parameters )64| . 
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these two curves is excellent as clearly illustrated in Fig.[T] Our approximate Sh{f) function can be easily calculated 
according to the following simple prescription. Define u ~ 'infr, where / is frequency and r = (5 x 10^ km)/c — 50/3 
sec is the light travel time down one of LISA's arms. For u < utrans = 0.25 we set 



Shif) 



.08 X 10-"^ 



5.52 X 10 



'41 



while for u > Wtia 



where 



1 / 2.88X10- 



-48 



5.52 X 10" 



/, 2/\a/1 2\ 4sin(u) cos(m) 
(l + cos2(u)) ( 3 - ^j +sm2(u) + 



(36) 



(37) 



(38) 



Note that LISA's characteristics are incorporated in the light travel time r and in the numerical constants in the above 
expressions, so a similar mission but with different noise characteristics would still be described by the functional form 
defined by eqns. ([5511551) . 

We have also included noise from a confused population of galactic white dwarf binaries following the prescription 
outlined in 13111. 



le-15 



le-16 



le-20 



le-21 



— LISA Average sensitivity, SNR=1 

— Analytic fit 




le-15 



0.0001 0.001 0.01 0.1 
frequency (Hz) 



- le-16 E 



le-17 



le-18 



le-19 



le-20 



le-21 




0.0001 0.001 0.01 0.1 
frequency (Hz) 



FIG. 1: Expected sensitivity curve ■s/KU) for LISA; black curve: numerical curve as generated in [63], red curve: analytic 
approximation used in this paper, see text for details). 



IV. KLUDGE WAVEFORMS AND FLUXES: RESULTS 



A. Comparison to Teukolsky-based geodesic waveforms 

1. Time-domain comparison 



In this Section we compare NK waveforms (using the overlap function defined in the previous Section) to a variety 
of TB waveforms from inclined-circular, equatorial-eccentric and some generic Kerr orbits. The TB waveforms are 
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treated here as "exact" (though this is not strictly the case as discussed in Section |n|. Our comparison rule is that 
the waveforms are generated using the same Kerr orbit with identical initial conditions. We label these orbits by the 
triad of orbital elements {p, e, i} defined in equations (I3|)-(I11). 

Figures [H [31 and S] show time-domain waveforms for a selection of orbits in the moderate and strong field. The 
NK waveforms plotted in these figures were generated by using the quadrupole-octupolc formula (|18p . The results for 
"plus" and "cross" polarizations are essentially the same, so in the figures we show only the "plus" polarization. For 
all of the figures, we have assumed an observation point located at $ = 0° and 8 = 45° or 90° (as indicated). 

Direct visual inspection of the waveforms gives some indication of how well the NK and TB waveforms match. In 
each of the cases illustrated, the kludge waveforms manage to capture the overall wave pattern, and for the orbits 
with p > 8M they almost exactly match the TB waveforms. For ultra- relativistic orbits (e.g., with p < 5M and 
1 — a/A/ ^ 1), the finer structure in the TB waveforms is clearly not reproduced by the quadrupole-octupole kludge 
waveforms. As we have discussed, these features are the imprints of higher multipole components in the radiation 
which are amplified by back-scattering, and thus are not expected to appear in the NK waveforms. Waveforms 
generated using the Press formula (not shown) do have some finer features due to the presence of higher multipole 
components, as we might have hoped, but they are not nearly as complicated as the structure of the TB waveforms. 

One also notices that for certain parts of the waveforms (e.g., Fig. there is a disagreement in the amplitude 
while the phase is accurately reproduced. This amplitude discrepancy is periodic, i.e., the points where the amplitude 
is poorly reproduced occur at regular intervals. This suggests that the kludge waveforms are missing some periodic 
components, as we might expect since they are truncated expansions in multipole moments. Indeed, the amplitude 
discrepancy is less pronounced in the Press waveform, which includes contributions at all multipoles. For the purposes 
of data analysis, it is important to have templates with a high overlap with the true signals, and to achieve that it 
is much more important that the waveform phase is reproduced rather than the waveform amplitude. The waveform 
phase is determined by the orbit generating the gravitational radiation. The fact that the kludge waveforms are based 
on true geodesic orbits is presumably the reason that we find, post facto, such impressively high overlaps with TB 
waveforms and especially for those from circular-inclined orbits. 

A comprehensive set of data for the overlaps between NK and TB waveforms is given in Tables |T1 |TT] and IIIII in the 
Appendix. These were computed using the overlap function described in Section [TTTl and assuming a central black hole 
mass of M = 10^ Mq. This was chosen since preliminary event rate estimates suggest the inspirals of ^ IOMq BHs 
into ^ 10^ Mq SMBHs will dominate the LISA detection rate [1^. These tables indicate that if the orbital periastron 
Tp ^ 5M, the overlap between TB waveforms and both the quadrupole-octupole and Press waveforms stays above 
~ 0.95. We also find that both these expressions have better performance than the pure quadrupole waveforms (fTBjl . 
but there is little difference between quadrupole-octupole and Press waveforms. The kludge and TB waveforms begin 
to deviate significantly for strong-field, ultra-relativistic orbits with Vp < AM, with the overlap dropping to ^ 50% 
for orbits that come very close to the horizon. Disappointingly, the Press waveforms do not seem to do much better 
than the quadrupole-octupole waveforms in this strong field regime, despite the inclusion of additional multipole 
components. The Press waveforms include only the "direct" contribution to the higher multipoles, i.e., the piece that 
arises from fast motion in fiat spacetime. In true EMRI GWs, the higher multipole contributions are significantly 
enhanced by "tail" effects, i.e., the backscattering of radiation from the background geometry. Since we are using 
the Press formula in flat space, we do not include this back-scattering enhancement. This is presumably why using 
the kludge Press waveforms does not significantly improve the overlap with TB waveforms. The Press waveforms 
do perform consistently better for circular orbits and weak-field eccentric orbits, but the difference between the two 
approaches is usually small (with the exception of a few cases which we discuss in more detail later). We conclude 
that the quadrupole-octupole waveform model is sufficient and there is not much gain from using the computationally 
more intensive Press formula. 

To summarize, we find that kludge waveforms are accurate — and very quick to generate — substitutes for TB 
waveforms for all orbits with periapse Vp ^ 5M . In Schwarzschild, the separatrix between bound and plunging orbits 
is at {rp)s = (6 -f 2e)/(l + e) M, so the kludge waveforms will be accurate throughout any inspiral with eccentricity 
at plunge e ^ 1/3. Computations of inspirals into Schwarzschild black holes To] indicate that, in many cases, the 
residual eccentricity at plunge will be small, so that kludge waveforms will be suitable for the majority of Schwarzschild 
inspirals. For retrograde orbits around Kerr black holes (90° < t < 180°), the minimal periapse is at even larger radii, 
so that even weaker restrictions must be imposed on the eccentricity at plunge for the kludge to be valid. In contrast, 
for prograde orbits (0° < 6 < 90°), an increased black hole spin allows stable orbits to exist much deeper in the 
strong field. For a = 0.9M, the separatrix for circular orbits is at rp « 2.3M, 2.6M, 3.7M for orbits with inclinations 
of 0°,30°,60° respectively; for orbits with a plunge eccentricity of e = 1/3, it is at rp « 2.0M, 2.2M, 3.1M for the 
same inclinations. As a ^ M , the separatrix of equatorial orbits decreases even more, asymptotically approaching 
(rp)^° = M [l^. Kludge waveforms are not particularly good in this regime, with overlaps ~ 50%; fortunately, this 
corresponds to a comparatively small region of parameter space. The evolution proceeds through this region very 
quickly, so we do not lose much signal-to-noise ratio by failing to match the waveforms in this region. 
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p=6.5M, e=0.5, a=0.5M, ej=90(deg) 




p=6.5M, e=0.5, a=0.5M, 9j=45(deg) I — TB waveform 
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p=3.5M, e=0.4, a=0.99M, e^=90(deg) I— NKwavefornil 




time (M) 

FIG. 2: Comparing TB and quadrupole-octupole kludge waveforms (black and red curves (bold black and grey curves in the 
b&w version) , respectively) for equatorial orbits and for an observer at a latitudinal position 9 — 45° or 90° . Orbital parameters 
are listed above each graph. The waveforms are scaled in units oi D/fi where D is the radial distance of the observation point 
from the source and n is the test-body's mass. The x-axis measures retarded time (in units of M) and we are showing the '+' 
polarization of the GW in each case. The overlaps between the NK and TB waveforms are 0.970, 0.987, 0.524 going from the 
top figure down. 

If the overlap between a given signal and the best-fit template in a search bank is less than 1 , this leads to a decrease 
in the maximum distance to which that signal can be detected, and a corresponding reduction in event rate. For the 
purposes of detection, overlaps as low as 50% might well be considered good enough, if the astrophysical event rate 
is sufficiently large [2C|. However, for the purposes of parameter estimation, much higher overlaps (> 95%) will be 
required in general. It is clear from the results in this paper that this is only partially achievable by the existing family 
of kludge waveforms. Nonetheless, these waveforms might be useful for LISA data analysis as search or detection 
templates over some (perhaps a large part) of the astrophysically relevant portion of the {a/M, p,e,t} parameter 
space. The waveforms may also provide sufficiently accurate estimation of the source parameters (in certain regions 
of parameter space) that they could be used as the first stage in a hierarchical search. The purpose of such a search 
would be to identify "interesting" regions of parameter space for follow up with more accurate waveforms. This is the 
main conclusion of this paper. 

However, the regions where kludge waveforms are good enough must be identified more carefully by comparison to 
accurate inspiral waveforms. As we have discussed earlier, the flat space emission formula used in the construction of 
the NK waveforms ignores all effects of scattering from the background curvature. These "tail" terms make a significant 
contribution to the waveform structure, and build up over the course of an inspiral. Although we have found good 
overlaps with geodesic waveforms and circular-inclined inspiral waveforms (see Section HV Bp . comparisons to inspirals 
of eccentric-inclined orbits are required to properly assess the importance of including the tail terms. Accurate, 
self-force waveforms for such orbits will not be available for a few years and only then will it be possible to firmly 
demarcate the regime of usefulness of the present, or further improved, NK waveforms. 




2. Frequency domain comparison 

To better understand the overlaps quoted in the tables we consider here the integrand of Eq. (|55)) for the value of 
t which maximizes the overlap: 

Hf)h*if) + i*{f)Hf)^^2.ft 



Shif) 



e""-^\ (39) 



In Figure[5]we plot this function for the generic orbit p = 12M, e = 0.1, i = 120° and a = 0.9M. For the signal i(/) 
we use the TB waveform and correlate it with templates h(f). As templates we use the TB waveform itself (black 
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r=5M, i=60(deg), a=0.99M, 9 =90(deg) 
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FIG. 3: Comparing TB and quadrupole-octupole Itludge waveforms (black and red curves (bold black and grey curves in the 
b&w version), respectively) for circular-inclined orbits and for an observer at a latitudinal position 9 — 90°. Orbital parameters 
are listed above each graph. The waveforms are scaled in units oi D/fi where D is the radial distance of the observation point 
from the source and fi is the test-body's mass. The x-axis measures retarded time (in units of AI) and we again show the plus 
polarization of the GW. The overlaps between the NK and TB waveforms are 0.882 for the top figure and 0.955 for the bottom 
figure. 



solid line), the Press waveform (squares), the quadrupole-octupole NK waveform (crosses) and the quadrupole NK 
waveform (circles). We have deliberately chosen a case for which including higher harmonics significantly (> 10%) 
improves the overlap, especially in the presence of white dwarf confusion noise (see Table |TTT]). 

One can see that the main contributions to the overlap come from several dominant harmonics. For a circular 
equatorial orbit the main harmonic would correspond to twice the orbital frequency, but in general the main harmonics 
depend on eccentricity 31]. Besides the harmonics coming from the azimuthal motion there are many additional 
components coming from beating between harmonics of the fundamental frequencies of the (j)-, 9- and r-motion: 
n^, ilgj^lr (see [g^ for a Fourier decomposition of the orbital motion). Depending on the mass M of the central BH, 
the dominant harmonics may lie in the most sensitive part of LISA's frequency range or they can be suppressed by 
confusion noise. In the latter case the higher harmonics, although smaller in amplitude, are effectively enhanced by 
the inverse power spectral density and therefore play an important role in SNR accumulation. Looking at the inlaid 
box in Figure [5] we can compare the quadrupole-octupole and Press waveforms at high frequency. Here, the Press 
waveform does perform better - e.g., around 3mHz the quadrupole-octupole NK waveform has failed to reproduce 
some harmonics in the TB signal, whereas the Press waveform does so pretty well. We will discuss the frequency 
representation of the signals further in the next subsection. 



B. Comparison to Teukolsky- based inspiral waveforms 

As described earlier, most of the existing TB waveforms describe GW emission from geodesic orbits, which is why we 
have focused on these for the comparisons above. However, TB waveforms have also been constructed for the inspiral 
of circular- inclined orbits [s^l- To illustrate the applicability of the kludge to the detection of realistic, inspiralling 
EMRI waveforms, we here compare a kludge inspiral waveform to a TB circular-inclined inspiral waveform generated 
using the method described in [37| . The kludge inspiral trajectory is constructed using the flux formulae described 
in [10], and then the particle trajectory is computed as discussed in Section fll C II We use the quadrupole-octupole 
formula (jlSp to generate a kludge waveform from this trajectory. 

A visual comparison is shown in Figure [HI There we plot both the inspiral trajectory r{t) and the gravitational 
waveform One can see an excellent agreement up to the last several cycles, where we would not trust either type 
of waveform anyway. We have used a large mass ratio (10~^) to compute the waveforms in order to speed up the 
orbital evolution. We should point out that the comparison is not improved if we use Press waveforms instead of 
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a=0.9M, p = 12M, e = 0.3, i = 140(deg), 6 = 60(deg) 
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FIG. 4: Comparing TB and quadrupole-octupole kludge waveforms (black and red curves (bold black and grey curves in the 
b&w version), respectively) for generic orbits. Orbital parameters are listed above each graph. The waveforms are scaled in 
units oi D/fi where D is the radial distance of the observation point from the source and /i is the test-body's mass. The x-axis 
measures retarded time (in units of A/). The overlaps between the NK and TB waveforms are 0.991 and 0.966 for the top and 
bottom figures respectively. 
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FIG. 5: Comparison of the integrand of the SNR in the frequency domain [h{f)s* (f)]/ Sh{f) where h{f) is determined from 
the quadrupole, quadrupole-octupole and Press waveforms, and s(/) is determined from the TB waveform. We have chosen 
orbital parameters p — 12M, e — 0.1, t = 120° and spin a — 0.9Af. One can see that Press waveform performs better than 
quadrupole-octupole waveform at high frequencies. 



quadrupole-octupole waveforms. 

We may quantify the agreement in the waveforms by again using the overlap function, which we plot as a function 
of truncation time in Figure [7] for three different masses of the central black hole. One can see that for the majority 
of the inspiral we obtain quite high overlap, but it decreases when we include the final stages of the inspiral. For the 
inspiral in Figure [HI the overlap (when M — 10^ A/q and in presence of the white dwarf background) up until radii 
of r = 8M, 6M, 5M, 4M, 3.5M is 0.985, 0.974, 0.92, 0.79, 0.755 respectively. The overlap over the whole inspiral is still 
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greater than 75%. An overlap of 75% corresponds to a factor of ~ 2 reduction in event rate. However, overlaps of 90% 
or 95% correspond to a reduction in event rate of only 1.4 or 1.17 respectively. A significant amount of signal to noise 
ratio accumulates over the early stage of inspiral, and so it may well be that kludge inspiral waveforms can be used 
to efficiently detect and track events until they reach r « 5 — 6M. Alternative techniques could then be employed for 
detection during the latter stages of inspiral. 

The mass of the central BH determines the frequency range through which the inspiral occurs. One could also 
evaluate the spectrum of the signal as a function of the dimensionless frequency (Mf) - the effect of changing M is 
then to shift the PSD Sh{f) right as M is increased or shift it left as M is decreased. In Figure [5] we again plot the 
overlap integrand ((39)l for AI — 10^ Mq. We also overplot the cumulative overlap 

^ Hf)h*{f)+S:*i.f)h{f) ,r 

and a scaled PSD. By shifting Sh{f) on the plot left/right one can see how the central mass affects the amount each 
harmonic contributes to the SNR/overlap. The higher harmonics are modeled less accurately by NK waveforms and 
this leads to an overall drop in the overlap for M — 10^A/q. However, the majority of the high frequency contribution 
to the waveforms comes from the end of the inspiral which we do not reproduce very accurately, rather than from 
higher harmonics emitted during the earlier, moderately relativistic part of the orbit. 

We should also mention that there is a problem with spectral leakage ^] (due to the large dynamic range) in the 
estimation of GW spectra. This is not a big problem on its own but for high M this leakage (at high frequencies) 
is amplified by the inverse power spectral density and leads to an erroneous result. In order to reduce this effect we 
have applied a window function 

w{t) = ^ (1 + tanh [A{t - nT)] + tanh [A{T - t - kT)]) 

to the template and signal for M = lO^AfQ and in the presence of confusion noise. The parameters k and A govern 
the steepness and cut-off point of the window; T is the waveform's duration. This function reduces the effect of the 
end of inspiral on the overlap and leads to an artificial increase in the overlap by a few percent. 

As mentioned above, a significant amount of signal to noise ratio accumulates over the early stages of the inspiral. 
To illustrate this we show in Figure [H] the accumulation of SNR as a function of time. To generate this plot, we have 
assumed the TB waveform is the signal and have used the NK waveform, truncated at different times, as the search 
template. One can see that we can recover up to 85% of the maximum SNR (i.e., the SNR if the template and signal 
were identical). Including the last few waveform cycles we see a drop in the SNR (^ 8%) due to the mismatch in the 
waveforms at the end of the inspiral. 

Employing Press waveforms helps to improve overlaps further, but for ultra-relativistic orbits the omitted tail 
contribution to the waveform becomes increasingly significant. 




C. Comparison to Teukolsky- based and PN fluxes 

It is possible to construct expressions for the flux of energy and angular momentum carried by the kludge GWs, 
and use these as estimates of the energy and angular momentum lost from the orbit in a true inspiral. The energy 
and angular momentum content of a TT GW field, h'i^rp, propagating in flat space at large distances from the source 
is determined from the Isaacson energy momentum tensor of the wave field p§i . [69| 

^GW-^{h^i^h^-). (40) 

Integration of this expression gives the energy and angular momentum loss rates due to GW emission (equations 
(4.13) and (4.22') of [t^) as 

^ = _i_ j (hi^^h^^Y''^^ (41) 
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FIG. 6: The orbital evolution, r{t), and h+ polarization of the GW signal, for an inspiralling quasi-circular inclined orbit with 
initial parameters p = lOM, e = 0, t = 45° and a = 0.9M, fi/M = 10"^, 6 = 90°. The black curve corresponds to the TB 
waveform and the red (grey in b&w version) line represents the quadrupole-octupole NK waveform. The inlaid box shows the 
part of the waveform where the TB and NK results start to deviate, which correlates with the deviation in the inspiral. 



M= 10 
6.52 5.96 5.15 3.19 



M=10^ 
6.52 5.96 5.15 3.19 



M= 10 
6.52 5.96 5.15 3.19 



0.95 



0.9 



0.85 - 



0.75 




I 

16800 



I 

18900 



0.9 




0.8 





h 

+ 


with WDs 


■ — ■ 


h 

+ 


no WD 


X— X 


h 

X 


with WD 




h 

X 


no WD 



16800 18900 
t/M 




16800 



18900 



FIG. 7: Overlap between TB and NK waveforms {h+ and hx) for an inspiralling (quasi)circular-inclined orbit as a function of 
truncation time. We show overlaps for different masses of the central BH: M = 10^, 10®, Mq and for sensitivity curves with 
and without white dwarf (WD) confusion noise. One can see that the overlap drops as we increase the truncation time and the 
mass M. On the top horizontal axis we show the orbital radius r/M corresponding to the truncation time. For A/ — W^Mq 
the lines with and without WD confusion lie on top of each other. 
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FIG. 8: The integrand of the overlap between TB and NK waveforms in the frequency domain (solid black line). Breaks in the 
curve correspond to a negative correlation. We also overplot the accumulative overlap (red solid line) and scaled Shif) as a 
dashed blue line. 
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FIG. 9: In this plot we present SNR obtained by applying the truncated NK waveform as a template to the full inspiral 
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normalised by the total SNR: SNR{s'^'^ , s'^^) = {s 



represented by TB waveform: SNR{ht^ ,0 j — ,^NKihNKMi2 

The primary mass was taken to be IQ^Mq and the sensitivity curved contained WD confusion noise. The total SNR for the 
source at 10 Gpc is approximately 66 (high SNR is due to large (10~^) mass ratio). 



Using the quadrupole-octupole formula to generate these expressions may be written in terms of the multipole 
moments of the source as 
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In these expressions J^j. = 



e.g.. 



Tka (2) 
^TT 



^jimS^'^ and bracketed numbers in superscripts denote the number of time derivatives, 
d^(/|.5')/dt^. If the kludge waveform is computed using the quadrupole formula only the leading 



order Ij^rp terms remain. If a Press waveform is used, there is an infinite sum of multipole components, and therefore 
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expressions pl1 - [42)) must be employed directly. 

The angle brackets in expressions (|41fB^ mean "average over several gravitational wavelengths" . This is achieved 
by taking a time average of the instantaneous flux expression inside the angle brackets, i.e., 



<X>-- 



Xdt 



(45) 



where T is an appropriate averaging time, ft is important to ensure that this averaging time is long enough. For 
circular equatorial orbits, the symmetries of the system ensure that the instantaneous and averaged fluxes are identical. 
For an eccentric equatorial orbit, there is only one periodicity in the flux — the period of the radial motion. Thus 
integrating over one radial period is sufficient to give an accurate average flux. For generic orbits, the averaging 
is more complicated since the periods of the motion in the radial and 6 directions are in general different and non- 
commensurate. The averaging time therefore has to be long enough to encompass several periods of both motions. We 
illustrate this in Fig.fTUlby plotting E as function of T for a generic orbit with parameters p = 4M, a = 0.99M, l = 60° 
and e = 0.1,0.4. The flux has been computed using the pure quadrupole waveforms ([16]). It is clear from Fig. [10] 
that the averaged flux converges over time, and this convergence is more rapid for the lower eccentricity orbit. In all 
subsequent flux calculations, we took T = lO^M (which corresponds to ~ 6 days for M — 10^ Mq). 

In Figure [TT] we show the angular distribution of the gravitational radiation from an EMRI orbit of given p and 
e, but for several different orbital inclinations, t. This picture is more or less typical of the majority of orbits. The 
variation of the energy flux with the sky position of the observer is indicative of beaming. This is strongest for 
equatorial orbits and reduces as the orbital inclination is increased. Equatorial orbits are restricted to a single plane, 
while inclined orbits wander through more of the spacetime. This wandering averages out the beaming over the sky. 
The larger the orbital inclination, the more of the phase space the body explores, the more averaging occurs and the 
more homogeneous the sky distribution of the energy flux. 

In Tables ITVI fc Fvl we show flux data for a series of orbits and black hole spins. We tabulate the flux computed via 
solution of the Teukolsky equation, and fluxes computed using quadrupole kludge waveforms and equation (|43p . We 
quote only the quadrupole results since the difference from using either the quadrupole-octupole or Press waveforms is 
quite small. As we mentioned in SectionJIT] Tanaka et.al. [2^ and more recently Gair et al. (38| computed approximate 
fluxes for Schwarzschild orbits based on the quadrupole waveform approximation (|16p . For a random choice of orbits 
we have found excellent agreement with their results. 

For further comparison, we also tabulate the "kludge" fluxes that have been used to construct inspirals and inspiral 
waveforms for scoping out LISA data analysis. The difference between these and our current results is a measure of 
the inconsistency in the kludged inspiral waveforms. In constructing an inspiral waveform, one set of expressions is 
used to evolve the orbit (given in [3, [l^ ) , and a different prescription (the one described here) is used to generate 
the waveform. This means that the energy carried by the kludge GWs is not equal to the energy lost by the orbit 
that is supposedly emitting these GWs. This leads to inaccuracies when the kludge waveforms are used to compute 
SNRs [131. In the original prescription for kludge inspirals [9], the phase space trajectories were generated using the 
leading order post-Newtonian results of Ryan [71i] . evaluated for relativistic orbital parameters. Ryan's expression for 
the energy flux is 



= (l-e^) 



2\3/2 



^^(')"]^(7j /2(e)cos6 



(46) 



where the e-dependent coefficients are, 
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In a more recent version of the kludge [l3| , this prescription for the evolution of the orbital parameters was improved. 
The new expression for E is complicated and we do not quote it here, but we also tabulate those results in Tables HV] 
& IV] under the heading "GG" . 

We note that in true black hole spacetimes, energy and angular momentum are lost into the horizon as well as to 
infinity. The Teukolsky equation provides a prediction for both the flux of radiation at infinity and the flux down the 
black hole horizon. When using Eqs. (|4ip ~ (|44p . we are only computing the energy and angular momentum carried to 
infinity in the kludge GWs. This should therefore only be compared to the Teukolsky flux at infinity. The discrepancy 
between these numbers is then a measure of the error associated with using the kludge waveforms to estimate signal 
to noise ratios for true inspirals. The TB fluxes quoted in the tables are accordingly just the infinity piece of the 
fluxes. When computing inspirals, which was the purpose for which the "GG" expressions were derived, the energy 
and angular momentum lost by the orbit should include the horizon piece of the flux. This should be born in mind 
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FIG. 10: Convergence of the averaged energy flux {J5) with respect to the averaging time T for a pair of equatorial-eccentric 
orbits. Note the increased convergence rate with decreasing eccentricity. 
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FIG. 11: Angular radiation pattern - energy radiated per unit sohd angle as a function of the colatitude of the observer, Qobs, 
for orbits with p = 5Af, e = 0.4, a = 0.9AI and a sequence of inclinations, as labelled. 

when considering the flux tables. However, the relative contribution of the horizon flux is only significant in the 
extremely strong field regime where both the "GG" expressions and the kludge fluxes cease to be valid. 

In general, we find that the fluxes computed from the kludge waveforms agree less well with Teukolsky results than 
do the "GG" fluxes. The kludge fluxes do improve over "GG" results for some retrograde orbits, and for some orbits 
with high eccentricity, but these are orbits for which the "GG" expressions are known to lack some terms. Once 
the "GG" expressions are further improved along the lines discussed in [l3|, it is likely that the "GG" fluxes will 
always be closer to TB results. Typically, the fluxes computed from the kludge waveforms underestimate the actual 
energy and angular momentum lost by the orbit, and also underestimate the "GG" fluxes. This suggests that inspiral 
waveforms that combine the two prescriptions will carry too little energy to infinity for the evolution seen in the orbit 
(or equivalently, the inspiral will proceed too quickly for the energy being carried by the GWs). In [38], a similar 
comparison was performed for inspirals into Schwarzschild black holes and they found that the kludge GWs contained 
too much energy, in contrast to this result. This difference arises primarily because, in that paper, inspirals were 
evolved using the "Ryan" fluxes, rather than the more up to date "GG" fluxes. The signal to noise ratio squared of a 
monochromatic source contributed when the frequency is in the range / — > /-|-d/ is proportional to Edf /{f^ f Sh{f))- 
The kludge waveforms get / largely correct (otherwise the high overlaps would not be maintained), but the error in 
E, arising from an amplitude discrepancy, will lead to an error in SNRs computed using kludge inspiral waveforms. 
Since E is in general an underestimate, the kludge SNRs will probably be underestimates of the true SNRs. The 
discrepancy in E is at most a few tens of percent. A discrepancy of 25% throughout the inspiral would lead to an 
error in the SNR of only 13%. However, the actual error will be much less, since for the majority of the inspiral the 
kludge E is much closer to the true value, and, as we saw earlier, most of the SNR comes from the early inspiral. 
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The kludge waveforms are thus accurate enough to be used for SNR calculations to estimate astrophysical event rates 
[20| . in which the astrophysical uncertainties far outweigh the waveform uncertainties. However, we must emphasize 
that amplitude discrepancies of this sort do not limit the applicability of the kludge waveforms for data analysis, 
provided the overlaps with true waveforms are high. If the kludge waveforms were used for parameter extraction, such 
discrepancies would lead to an error in the distance estimate only. In practice, we are likely to use kludge waveforms 
for detection rather than source characterization, and for that the waveform amplitude is irrelevant. 

Table IVl tabulates data for "parabolic" orbits, e = 1. The GW emission from such orbits is important for the 
capture problem, i.e., the mechanism by which compact objects are put onto EMRI orbits. Gravitational radiation 
from such orbits has been studied in the past using the Teukolsky formalism [40|,[6H, but tabulated data has only 
appeared in a recent time-domain analysis by Martel [l^]. Those results are claimed to be accurate at the level of 
< 1%, and therefore can be treated as exact for our purposes. Table IVl contains the relevant data together with our 
kludge results. For a periastron fp ~ 5M the error in using the kludge fluxes is ~ 30% and rapidly decays (grows) 
as we move to larger (smaller) periapse. In this Table we have also included some data for non-equatorial parabolic 
orbits, but for this class of orbits there are as yet no available TB results. Our experience with bound generic orbits 
suggests it is likely these data will have accuracies of the same order as the equatorial orbits with a similar periastron. 
In Fig. [T^l we plot AE as a function of Vp and (left and right panels, respectively) for parabolic orbits with 
a/M = 0, 0.99 and l — 0°, 60°, 120°, 180°. All curves are computed using the quadrupole waveform flux p6| and 
are terminated at the point where the orbit plunges. Fig. [T^ has astrophysical significance as it displays (within the 
accuracy of the present calculation) the amount of energy lost in a single parabolic (or e ~ 1) encounter with a SMBH, 
for a variety of inclinations and for a range of periastra that are appropriate to the capture problem (sj, [ill- This 
energy loss data can be used to estimate the capture rate of compact objects by a single massive black hole. 

We note that it is clear from Table |V] that the "GG" results perform quite poorly for parabolic orbits. This is 
because the "GG" results are built on a small eccentricity expansion, which is no longer valid when e = 1. At present 
the kludge fluxes are closer to TB results than the "GG" fluxes where comparisons can be made, and most likely 
are better approximations for the other orbits as well. The nerformace of the "GG" fluxes for highly eccentric orbits 
is now being improved, ; oximate TB fluxes 
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FIG. 12: Total radiated energy for Kerr parabolic orbits; as function of periastron. 



V. INCLUSION OF CONSERVATIVE CORRECTIONS 



In the preceding sections, we have shown that kludge waveforms do very well at reproducing the GW emission from 
geodesic orbits, and, if coupled to an accurate prescription for the radiation fluxes, can very accurately reproduce 
the GW emission from inspirals. However, the radiation fluxes are not the whole story, as they represent only the 
dissipative piece of the self-force. In true inspirals, the self-force also has a "conservative" component. This does 
not lead to evolution of the orbital parameters and hence cannot be determined from the flux of energy, angular 
momentum and Carter constant, but does lead to a change in the phasing of the gravitational waveform. There has 
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been some debate recently about the importance of including the conservative pieces of the self- force when constructing 
templates for EMRI source detection [71, [tI] . It is quite possible that these corrections will be needed for detection, 
at least for some range of source parameters; these terms will certainly be needed for parameter determination. We 
aim to demonstrate here how easily conservative terms can be included in the kludge waveform prescription. 

The conservative part of the self-force will include an oscillatory component, which averages to zero, and a secular 
piece that leads to accumulation of a phase error over time. The net effect of the secular piece is that, for a given set 
of orbital parameters {p, e and t), the average frequencies of the motion in the r, 9 and (j) directions will differ from 
the expected values for a Kerr geodesic. The conservative self-force might also affect the harmonic structure of the 
GWs, but the dominant contribution will come from orbital dephasing. In the spirit of the kludge, we will now try 
to account for this part of the effect. 

The effect of the conservative self-force can be characterized by the changes in the r, 9 and (j) frequencies, averaged 
over many orbits. We denote these frequency changes by SQ^, and (5f20. Conservative effects can then be included 
in the kludge orbital determination simply by amending the evolution equations for the angular variables ip, x s^nd (p 
given earlier (Eqs. ([8]), (l9|), and (fTc]) respectively): 

(^)^(^)..— (^)^(^)./--- 

The subscript "geo" denotes expressions for the phase derivatives in geodesic orbits given earlier. In general, the 
frequency shifts will depend on the orbital parameters, and for an inspiralling source we assume that these are 
evaluated for the instantaneous orbital parameters, as discussed in Section [II C 11 

To evaluate these frequency shifts in the framework of black hole perturbation theory requires self-force calculations 
[2^ which are not yet fully developed, but should be on the timescale of a few years. Once the self-force corrections 
are known, it should be possible to use them to compute the effective frequency shifts in the framework we use here, 
i.e., objects moving on geodesies of the Kerr spacetime, the parameters of which evolve with time. In the meantime, 
conservative corrections are already known in the post- Newtonian framework up to 3.5PN order pTSj. These can be 
used to compute leading order conservative corrections to include in the kludge. We outline here a general method 
for such a calculation, illustrated by the simplest case of a IPN conservative correction to circular orbits in the 
Schwarzschild spacetime. The generalisation to 3.5PN corrections for arbitrary orbits in Kerr will be presented in a 
future paper. 

The greatest difhculty in such a computation is to identify coordinates between the post-Newtonian and perturbative 
formalism. The former calculations tend to be carried out in harmonic coordinates while the natural coordinate 
system to use for the latter are black-hole centered coordinates, like the Boyer-Lindquist coordinate system. The best 
approach is to use asymptotic observables, specifically the orbital frequency, perihelion precession frequency, orbital 
plane precession frequency and their first time derivatives, to identify coordinates between the two formalisms and find 
the missing conservative pieces. The power of this approach is that it corrects the kludge in a physically meaningful 
way; however, this comes at the cost of an effective 0{rj) redefinition of the orbital parameters, where 77 = fi/AI. For 
extreme mass ratios, rj < 10~^, this is less than the expected error in parameter determination [3l| . 

The post-Newtonian model has two particles moving on orbits (expanded to the stated post-Newtonian order) and 
evolving via post-Newtonian radiation reaction expressions. These orbits include conservative corrections to geodesic 
motion. For the case of circular orbits in Schwarzschild we find at IPN and with mass ratio corrections at linear order 
only that the orbital frequency and its time derivative are given by (76j 

dt \Rj M \ 2 \Rj \2 i 'J J 

df7 _ 96 77 f^'^\^ fi 3 /M\ /2591 13571 

'dt ~ Tap \r) \ ^2^^ yRj ysm ^ 672 ^ 

In the above, we use M to denote the mass of the primary (not the total mass as is conventional in the post-Newtonian 
approach), and R to denote the post-Newtonian orbital semi-major axis. In the kludge, the orbit is a geodesic of 
the Kerr spacetime, and the orbit is evolved according to the prescription in [l3|. For circular equatorial orbits in 
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Schwarzschild we have (at IPN order) 



1 { M\ , , 

O = - — (51) 



r 



At 5 M \ r J \ 336 \ r J J ^ ^ 

an 96 7] fM\^ f 743 fM\ ^/ 3\\ , , 

— l-TT^ — +0 (53) 



dt 5 AP \ r J \ 336 V ^ . 
Here, r denotes the Schwarzschild radial coordinate of the orbit. If we write 

r = i?. (^1 + (^^^ ba + ■■■+!] CQ + rj (^^^ ci + • • • ^ (54) 

and substitute into expression (|5ip or (|53p . keeping corrections at IPN order and linear in the mass ratio, we can 
compare to equations or ([50]) and solve for bo, cq and Ci. The comparison between ([ST]) and gives a different 
result for cq and ci than the comparison between (|53p and (j50p . since the kludge prescription ignores conservative 
effects. We therefore write 

and solve simultaneously for 5o, Cq, Ci, do di [bearing in mind that (j52p not (|53p is the fundamental quantity to 
examine for computing inspiral in the kludge]. We find that the required parameters are 

K ^ 1 845 1 1975 

On = 1, Cn = , Ci = , Ctn = — , "1 = ■ 50 

" ' "4' 448' 8' 896 ^ ^ 

We note that the correction to Q scales as R~^^^, compared to the scaling of tl, which is R^^^^^. In the weak-field, 
i? — !■ 00, the conservative correction will be much more significant than the radiative correction. This is in keeping 
with what Pound et al. [74| found for a particle moving under an electromagnetic self-force. 

Extension to generic orbits is straightforward. In that case, we have two extra observables, the perihelion precession 
frequency and the rate of precession of the orbital plane, and two extra coordinate dependent orbital parameters, the 
eccentricity and orbital inclination. Identifying the rate of change of the orbital frequency, the perihelion precession 
frequency and the orbital plane precession frequency between the post-Newtonian and kludge approaches gives us a 
relation between the two sets of coordinates. Comparing the values of the frequencies then tells us which Jfi^, 5 fig 
and Sflr must be added into the kludge. In all of this we have chosen implicitly to identify the two masses, M and 
fi, and the spin, a, between the post-Newtonian and kludge models. To our minds, it makes more sense to hold these 
parameters fixed between the two models, at the cost of modifying the definitions of the orbital elements p, e, l. 

Once self-force data is available, it will be possible to obtain a consistent solution for the conservative correction 
in Boyer-Lindquist like coordinates by using fits to the self- force results. In using equation ([5^ . we have implicitly 
assumed the relationship between energy/angular momentum and the radial coordinate is unchanged by the self-force. 
However, in self- force calculations, the energy and angular momentum of circular orbits at a given radius do change 
as a result of conservative effects. This is merely a manifestation of the 0{r]) redefinition of r implicit in this approach 
that we mentioned earlier. When comparing to self-force calculations, it will be preferable to include this redefinition 
of energy as well as the conservative corrections to Q so that r maintains it's meaning. Ultimately, both approaches 
are equivalent. Another way to obtain cIq and di for equation (j55p is by expanding dQ/dt as a function of Q in 
both the PN and kludge approaches. Matching various orders in this expansion gives do and di directly without 
having to simultaneously solve for the coordinate transformation. The function dfl/dt{n) is a GW observable, so 
this approach ensures that the kludge will have the correct leading order form for this observable. However, this is a 
series truncated at a given PN order, so it will be important to assess how significant the omitted terms can be. It is 
likely that using the present kludge, augmented with conservative corrections up to 3.5PN order computed as outlined 
above, will do well enough at reproducing the phasing of true inspiral waveforms. As discussed in |74l |. it is likely that 
the conservative terms will be most important in the weak-field regime, where the post-Newtonian results are valid. 
In the regime where the post-Newtonian corrections cease to be valid, the contribution of conservative corrections to 
the phase evolution may be much less critical. This will be investigated in future work. If the PN expansion is not 
sufficient, the kludge can be augmented using fits to the results of self-force calculations as described. 
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VI. CONCLUDING DISCUSSION 

In this paper we have provided a simple, "easy-to-use" , prescription for approximating the gravitational waveforms 
generated by test-bodies inspiralling into Kerr black holes. These "kludge" waveforms are constructed by combining 
familiar flat-spacetime wave equations with a true geodesic trajectory in the Kerr spacetime. Despite its formal 
inconsistency this hybrid approximation results in remarkably accurate waveforms, which we have established by 
comparison to rigorous Teukolsky-based perturbative waveforms. We find an impressive overlap between the kludge 
and TB waveforms for particles on geodesic orbits (i.e., ignoring radiation-reaction). This overlap is > 95% over 
a significant portion of the relevant orbital parameter space. Significant degradation (overlap reduced to ~ 50%) 
occurs for strong field orbits around rapidly spinning black holes. For such cases, the contribution of the radiation 
backscattered from the background spacetime is sizable; this effect is not included in the present formulation of the 
kludge waveforms. As a rough (but reliable) rule of thumb, we have found that kludge waveforms work well for Kerr 
orbits with periaspe distance Vp ^ 5Af, irrespective of the black hole spin. 

We have experimented with three different types of kludge waveforms generated using three different solutions of 
the flat-spacetime gravitational wave-equation: quadrupole, quadrupole-octupole and "Press" waveforms [s^. The 
latter choice includes contributions from all multipole moments of the orbiting body. We have concluded that for 
all practical purposes the quadrupole-octupole waveforms are the optimal choice, nearly as accurate as the Press 
waveforms while much quicker and easier to generate. 

Within the adiabatic approximation, we have applied our method to calculate full inspiral kludge waveforms taking 
into account orbital evolution. This realistic scenario requires an additional "kludge" for describing the orbital 
evolution itself. Our earlier work [!i|,[l3| provides such a scheme: exact geodesic orbital dynamics coupled with 
approximate PN-based expressions for the fluxes. The resulting waveforms still have very high overlaps with available 
Teukolsky-based inspirals for circular-inclined orbits. The overlap is > 75% even when the inspiral terminates in 
the strong field region. This result indicates the kludge will be a very useful tool for generating inspiral waveforms. 
However, as yet Teukolsky-based inspirals are available only for circular-inclined orbits. Only when comparisons 
for generic inspiral orbits are possible will we be able to firmly identify the range of applicability of these kludge 
waveforms. 

We have also used the kludge waveforms to estimate the fluxes of energy and angular momentum carried away in 
GWs from geodesic orbits. We have found that such kludge fluxes in general are not as accurate for evolving inspirals 
as the post- Newtonian based expressions described in ^] . The kludge fluxes do provide us with an estimate of the 
error in using kludge waveforms for SNR calculations, and our results suggest the kludge waveforms will tend to 
underestimate the true SNRs if anything. 

The area where these kludge waveforms will find most use is in the development of EMRI data analysis for LISA. 
The combination of accuracy and simplicity of generation has already made these waveforms invaluable tools for the 
study of data analysis issues. It seems quite plausible that the waveforms may also play a role in the final search of 
the LISA data. One use could be for estimation of the waveform parameters as the first stage in a hierarchical search. 
The high faithfulness of these waveforms suggest that they may be able to set fairly tight bounds on the parameters 
of the emitting system. This will be extremely useful input for the second stage of the search where the system 
parameters will be refined using more accurate waveforms. However, as we discussed before, the kludge waveforms 
at present do not include some important physical features that we expect in true inspirals. Only when accurate 
Teukolsky or self-force based waveforms are available for generic inspiral orbits will we be able to fully quantify the 
range of validity and level of accuracy of the kludge waveforms. 

Time-frequency searches [tJ [zl; Um could also be used to detect EMRIs in the first stage of a hierarchical search. 
These can typically detect EMRI events at about half the distance of the semi-coherent search [23| , without the use 
of templates. A template with overlap of ^ 50% with the true signal can detect that signal at about half the distance 
of the correct template. The NK waveforms can easily achieve overlaps of ~ 50% through most of parameter space, 
particularly given the 2 — 3 week integration times needed for the semi-coherent approach. It is clear therefore that 
a first stage semi-coherent search using kludge waveforms would have a greater reach than time-frequency methods, 
and would put tighter constraints on the parameters. However, this would come at a much greater computational 
cost. The final LISA search will undoubtedly employ both methods at various stages in the analysis. 

A further possible application of kludge waveforms is to the study of non-Kerr EMRIs. It is hoped that LISA 
observations will allow "spacetime-mapping" of black holes [iO|j and thereby test the no-hair theorem. To 

carry out such tests quantitatively will require waveform templates, which incorporate the deviation from Kerr in 
the spacetime structure as a set of suitable parameters (e.g., multipole moments [l6[). The development of rigorous 
non-Kerr EMRI waveforms is a very difficult task, since generic spacetimes lack some symmetries that allowed the 
formulation of the Teukolsky framework for the Kerr spacetime. One way to make progress is to construct kludged 
waveforms in non-Kerr spacetimes along the lines outlined in this paper. Our results for the Kerr spacetime suggest 
that such kludged waveforms may well be sufficiently accurate for qualitative, and even quantitative studies of non- 
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Kerr EMRIs. All that is required is the integration of geodesic equations in the non-Kerr spacetime and thus kludge 
waveforms provide a computationally quick and easy tool to study "bumpy" and "quasi-Kerr" spacetime mapping 

The kludged waveforms/inspirals presented in this and a companion paper 10] could be considered as "second 
generation" , an improvement of the original simple version of quadrupole waveforms ^30j] and inspirals [9] . Certainly, 
there is space for further improvement. The inclusion of conservative self- force effects on the inspiral is an obvious next 
step, and we have already discussed here how this can be achieved by calculating the relevant orbital frequency shifts 
for circular Schwarzschild orbits. We can also further improve waveform generation by including the back-scattering 
effect arising from the propagation of the GWs in a curved spacetime. 

In summary, given the present level of performance of our kludge waveforms and inspirals and their prospects of 
improvement, we feel that they should remain valuable tools for LISA source modelling for the coming years. 
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TABLE L Numerical data for overlaps between TB and kludge waveforms - equatorial-eccentric Kerr orbits. Data is not shown 
with WD confusion noise for the ten orbits in the very strong field regime. In this regime, none of the kludge waveforms 
reproduce the TB waveforms very well, and this is compounded when the dominant harmonics are suppressed by white dwarf 
confusion noise. The overlaps are uniformly poor and uninformative, so we do not include them. 



p/M " 




(deg) 


a/M'' 


e= (deg) 


Quad 


overlap (+)^ 
Quad-Oct 


Press 


overlap with WD (+) 
Quad Quad-Oct Press 


duration(M)'' 


1.7 


0.1 





0.99 


90 


0.84 


0.772 


0.741 






2000 


1.7 


0.3 





0.99 


90 


0.76 


0.557 


0.500 






700 


1.9 


0.5 





0.99 


90 


0.544 


0.570 


0.547 






700 


1.9 


0.5 





0.99 


90 


0.523 


0.484 


0.445 






2000 


2.11 


0.7 





0.99 


45 


0.562 


0.566 


0.562 






700 


2.2 


0.7 





0.99 


90 


0.526 


0.496 


0.458 






700 


2.5 


0.1 





0.99 


90 


0.906 


0.853 


0.827 






2000 


2.5 


0.5 





0.99 


90 


0.671 


0.665 


0.651 






700 


3.5 


0.4 





0.99 


90 


0.588 


0.524 


0.507 






2000 


3.5 


0.4 





0.99 


45 


0.624 


0.598 


0.593 






5000 


5.1 


0.5 





0.5 


90 


0.856 


0.962 


0.967 


U.oOD 


U.yDl U.yDo 


700 


5.5 


0.5 





0.5 


90 


0.864 


0.964 


0.973 


0.862 


0.962 0.973 


2000 


6.0 


0.4 





0.5 


90 


0.871 


0.970 


0.980 


0.864 


0.967 0.979 


2000 


6.0 


0.5 





0.5 


90 


0.858 


0.966 


0.974 


0.855 


0.963 0.973 


2000 


6.5 


0.5 





0.5 


90 


0.870 


0.970 


0.979 


0.864 


0.968 0.978 


2000 


6.5 


0.5 





0.5 


45 


0.937 


0.987 


0.990 


0.932 


0.986 0.990 


8000 


10.0 


0.3 


180 


0.99 


90 


0.864 


0.961 


0.966 


0.806 


0.943 0.954 


8000 


10.0 


0.3 


180 


0.99 


45 


0.922 


0.971 


0.969 


0.883 


0.957 0.956 


8000 


10.4 


0.5 


180 


0.99 





0.998 


0.998 


0.999 


0.997 


0.999 0.998 


2000 


10.5 


0.5 


180 


0.99 


90 


0.878 


0.975 


0.982 


0.856 


0.968 0.978 


2000 


15.0 


0.4 





0.5 


90 


0.824 


0.963 


0.968 


0.60 


0.878 0.881 


8000 


15.0 


0.4 





0.99 


90 


0.824 


0.961 


0.963 


0.603 


0.874 0.860 


8000 



"semi-latus rectum 
''eccentricity 
'^inclination angle 
"^spin 

'^observation point, (f) = always 

■'overlap between "+" polarization of TB waveform with quadrupole ("Quad"), quadrupole-octupole ("Quad-Oct") and Press ("Press") 
kludge waveforms 
f over lap between "x" polarizations 
''waveform duration 
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TABLE II: Numerical data for overlaps between TB and kludge waveforms - inclined- circular Kerr orbits. 



p/M 


e 


i (deg) 


a/M 


e (deg) 




overlap (+) 






overlap(x) 




over 


ap with WD (x) 


duration(M) 












Quad 


Quad-Oct 


Press 


Quad 


Quad-Oct 


Press 


Quad 


Quad-Oct 


Press 




5.0 





30 


0.5 


45 


0.944 


0.984 


0.990 


0.946 


0.984 


0.990 


0.945 


0.984 


0.990 


3000 


5.0 





30 


0.5 


90 


0.899 


0.974 


0.984 


0.891 


0.971 


0.984 


0.89 


0.969 


0.983 


3000 


5.0 





30 


0.99 


45 


0.929 


0.969 


0.975 


0.93 


0.969 


0.975 


0.927 


0.967 


0.974 


3000 


5.0 





30 


0.99 


90 


0.904 


0.958 


0.964 


0.9 


0.957 


0.966 


0.895 


0.954 


0.965 


3000 


5.05 





60 


0.5 


45 


0.924 


0.967 


0.973 


0.922 


0.969 


0.977 


0.922 


0.968 


0.976 


3000 


5.05 





60 


0.5 


90 


0.91 


0.955 


0.961 


0.911 


0.963 


0.973 


0.907 


0.961 


0.972 


3000 


5.0 





60 


0.99 


45 


0.857 


0.912 


0.917 


0.86 


0.917 


0.925 


0.857 


0.913 


0.921 


3000 


5.0 





60 


0.99 


90 


0.854 


0.882 


0.888 


0.87 


0.912 


0.925 


0.862 


0.907 


0.922 


3000 


10.0 





30 


0.5 


45 


0.93 


0.989 


0.995 


0.936 


0.99 


0.994 


0.751 


0.957 


0.975 


8000 


10.0 





30 


0.5 


90 


0.89 


0.981 


0.990 


0.901 


0.981 


0.990 


0.63 


0.93 


0.969 


7000 


10.0 





30 


0.99 


45 


0.92 


0.98 


0.986 


0.93 


0.98 


0.986 


0.674 


0.92 


0.946 


8000 


10.0 





30 


0.99 


90 


0.884 


0.974 


0.982 


0.915 


0.975 


0.983 


0.371 


0.892 


0.963 


7000 


10.0 





45 


0.7 


45 


0.932 


0.980 


0.981 


0.936 


0.983 


0.987 


0.75 


0.935 


0.953 


8000 


10.0 





45 


0.7 


90 


0.922 


0.972 


0.972 


0.92 


0.98 


0.987 


0.68 


0.929 


0.969 


8000 


10.0 





60 


0.5 


45 


0.914 


0.981 


0.987 


0.912 


0.982 


0.990 


0.691 


0.937 


0.969 


8000 


10.0 





60 


0.5 


90 


0.945 


0.98 


0.981 


0.932 


0.981 


0.986 


0.732 


0.925 


0.948 


7000 


10.0 





60 


0.99 


45 


0.873 


0.953 


0.958 


0.875 


0.957 


0.965 


0.611 


0.875 


0.906 


8000 


10.0 





60 


0.99 


90 


0.923 


0.941 


0.937 


0.915 


0.959 


0.965 


0.675 


0.86 


0.893 


7000 


20.0 





30 


0.5 


45 


0.934 


0.987 


0.982 


0.938 


0.989 


0.989 


0.957 


0.987 


0.902 


50000 


20.0 





30 


0.5 


90 


0.889 


0.970 


0.922 


0.893 


0.974 


0.962 


0.92 


0.984 


0.900 


30000 


20.0 





30 


0.99 


45 


0.93 


0.98 


0.951 


0.94 


0.988 


0.992 


0.951 


0.987 


0.977 


50000 


20.0 





30 


0.99 


90 
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0.968 


0.907 


0.914 


0.979 


0.988 


0.933 


0.979 


0.951 


30000 


20.0 





60 
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45 


0.92 
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0.966 


0.916 


0.981 


0.995 


0.952 


0.958 


0.990 


50000 


20.0 





60 


0.5 


90 


0.95 


0.973 


0.958 


0.934 


0.976 


0.972 


0.95 


0.967 


0.957 


30000 


20.0 





60 


0.99 


45 


0.895 


0.963 


0.962 


0.892 


0.967 


0.973 


0.934 


0.951 


0.896 


50000 


20.0 





60 


0.99 


90 


0.954 


0.96 


0.938 


0.937 


0.972 


0.978 


0.963 


0.946 


0.943 


30000 
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TABLE III: Numerical data for overlaps between TB and kludge waveforms - generic Kerr orbits. 



p/M e b (deg) a/M Q (deg) 



overlap(+) 
Quad Quad-Oct 



overlap(x) 
Press Quad Quad-Oct 



Press 



overlap with WD ( x ) 
Quad Quad-Oct Press 



duration(M) 



6.0 0.1 
6.0 0.1 
6.0 0.1 



20.1364 0.9 
20.1364 0.9 
20.1364 0.9 



90 
60 
30 



0.912 
0.935 
0.972 



0.982 
0.988 
0.995 



0.993 
0.996 
0.998 



0.91 
0.939 
0.972 



0.984 
0.99 
0.995 



0.996 
0.997 
0.998 



0.894 
0.929 
0.968 



0.98 
0.989 
0.995 



0.995 
0.997 
0.998 



6.0 0.5 
6.0 0.5 
6.0 0.5 



20.103 0.9 
20.103 0.9 
20.103 0.9 



90 
60 
30 



0.89 
0.916 
0.958 



0.967 
0.975 
0.985 



0.973 
0.980 
0.985 



0.89 
0.912 
0.959 



0.972 
0.978 
0.985 



0.980 
0.982 
0.986 



0.875 
0.913 
0.955 



0.968 
0.976 
0.984 



0.979 
0.981 
0.985 



6.0 
6.0 
6.0 



0.1 
0.1 
0.1 



60.1461 0.9 
60.1461 0.9 
60.1461 0.9 



90 
60 
30 



0.951 
0.958 
0.943 



0.989 
0.989 
0.986 



0.993 
0.994 
0.996 



0.941 
0.95 
0.942 



0.988 
0.988 
0.987 



0.995 
0.995 
0.996 



0.933 
0.946 
0.941 



0.987 
0.987 
0.986 



0.995 
0.996 
0.996 



6.0 0.5 
6.0 0.5 
6.0 0.5 



60.1108 0.9 
60.1108 0.9 
60.1108 0.9 



90 
60 
30 



0.934 
0.939 
0.917 



0.975 
0.974 
0.971 



0.980 
0.982 
0.978 



0.919 
0.927 
0.916 



0.974 
0.972 
0.971 



0.982 
0.984 
0.979 



0.912 
0.924 
0.915 



0.971 
0.97 
0.971 



0.981 
0.984 
0.979 



6.0 0.7 
6.0 0.7 
6.0 0.7 



60.0755 0.9 
60.0755 0.9 
60.0755 0.9 



90 
60 
30 



0.926 
0.928 
0.897 



0.966 
0.967 
0.963 



0.970 
0.971 
0.966 



0.911 
0.919 
0.897 



0.967 
0.966 
0.963 



0.974 
0.974 
0.967 



0.906 
0.916 
0.897 



0.966 
0.966 
0.963 



0.975 
0.974 
0.968 



12.0 0.1 
12.0 0.1 
12.0 0.1 



119.9586 0.9 
119.9586 0.9 
119.9586 0.9 



90 
60 
30 



0.916 
0.926 
0.9 



0.987 
0.988 
0.983 



0.994 
0.991 
0.994 



0.911 
0.92 
0.899 



0.988 
0.989 
0.986 



0.998 
0.998 
0.994 



0.627 
0.634 
0.577 



0.922 
0.926 
0.91 



0.985 
0.985 
0.956 



12.0 0.5 
12.0 0.5 
12.0 0.5 



119.9686 0.9 
119.9686 0.9 
119.9686 0.9 



90 
60 
30 



0.938 
0.942 
0.924 



0.988 
0.988 
0.984 



0.992 
0.993 
0.993 



0.934 
0.936 
0.924 



0.988 
0.987 
0.986 



0.995 
0.994 
0.994 



0.889 
0.889 
0.862 



0.979 
0.977 
0.971 



0.993 
0.991 
0.988 



12.0 0.7 
12.0 0.7 
12.0 0.7 



119.9786 0.9 
119.9786 0.9 
119.9786 0.9 



90 
60 
30 



0.935 
0.940 
0.926 



0.986 
0.987 
0.983 



0.992 
0.992 
0.992 



0.933 
0.934 
0.924 



0.987 
0.985 
0.983 



0.993 
0.992 
0.990 



0.904 
0.904 
0.884 



0.981 
0.978 
0.973 



0.991 
0.990 
0.986 



12.0 0.3 
12.0 0.3 
12.0 0.3 



139.9597 0.9 
139.9597 0.9 
139.9597 0.9 



90 
60 
30 



0.9 
0.938 
0.948 



0.985 
0.991 
0.992 



0.992 
0.994 
0.994 



0.891 
0.931 
0.949 



0.984 
0.99 
0.993 



0.996 
0.997 
0.998 



0.751 
0.827 
0.859 



0.955 
0.972 
0.979 



0.991 
0.993 
0.994 



12.0 0.1 
12.0 0.1 
12.0 0.1 



159.9728 0.9 
159.9728 0.9 
159.9728 0.9 



90 
60 
30 



0.836 
0.881 
0.947 



0.973 
0.984 
0.994 



0.992 
0.995 
0.997 



0.838 
0.889 
0.947 



0.974 
0.987 
0.995 



0.993 
0.998 
0.997 



0.494 
0.549 
0.726 



0.887 
0.926 
0.96 
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TABLE IV: A comparison of the approximate fluxes computed from the kludge quadrupole waveforms (|16|) described in this 
paper to accurate Teukolsky based fluxes for circular (top section), equatorial (middle section) and generic (bottom section) 
orbits. The energy fluxes are expressed in units of IJ? /M, and the angular momentum fluxes in units of /^^/M^. The subscripts 
"TB", "GG" and "quad" denote fluxes computed from Teukolsky based waveforms, from the kludge flux expressions in [iy| 
and the approximate quadrupole fluxes respectively. 
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TABLE V: As in Table HV] but for parabolic orbits. 
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